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I.  INTRODUCTION 

The  U.S.  Navy  has  been  greatly  interested  in  determining  the  effects  of 
underwater  explosions  on  the  structural  dynamics  of  ships  and  submarines  since  well 
before  World  War  II.  Until  recently,  ship  shock  trials  and  experimental  tests  on  scale 
models  had  been  the  preferred  means  to  obtain  the  required  information.  However,  due  to 
the  high  cost  and  environmental  problems  associated  with  underwater  explosion  tests  and 
the  tremendous  increases  in  computer  processing  capability,  computer  simulations  have 
become  a  widely  accepted  alternative  to  the  experimental  tests  in  determining  the 
dynamic  response  and  failure  modes  of  submerged  and  partially  submerged  structures. 
These  simulations  can  analyze  scenarios  which  are  far  too  dangerous  to  attempt  with  ship 
shock  trials.  Furthermore,  they  provide  valuable  insight  into  the  design  and  shock 
qualification  process  resulting  in  more  combat  effective,  survivable  and  cost  effective 
ships  and  submarines. 

Composite  materials  are  beginning  to  see  widespread  use  as  major  structural 
components  of  military  vehicles.  They  are  typically  used  to  improve  the  stiffiiess-to- 
weight  ratio  or  the  strength-to-weight  ratio  of  a  structural  member.  The  mechanical 
properties  can  be  tailored  based  on  the  design  requirement  in  a  manner  unavailable  with 
conventional  metals  or  alloys.  While  manufacturing  difficulties  and  cost  have  limited 
their  use  in  current  ship  and  submarine  hulls,  they  have  great  potential  to  significantly 
improve  the  performance  of  future  classes  of  ships  and  submarines. 

The  objective  of  this  study  is  to  develop  a  general  purpose  finite  element  code  for 
the  structural  analysis  of  composite  materials  subjected  to  either  underwater  shock  or 
impact  loading  for  use  on  the  personal  computer  (PC)  platform.  Linear  elastic  behavior  of 
the  material  is  assumed  for  this  study.  At  a  later  time,  nonlinear  effects  can  be  added  such 
as  the  damage  progression  in  composite  materials  for  example. 

A  brief  overview  of  the  finite  element  formulation  used  in  this  code  is  presented 
first  followed  by  the  development  and  solution  of  the  fluid-structure  interaction 
equations.  Several  code  verification  problems  are  then  demonstrated  followed  by  an 
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application  of  the  code  to  a  sandwich  composite  impact  problem  and  a  study  of  the 
effects  of  layer  smearing  on  the  fluid-structure  interaction. 
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II.  FINITE  ELEMENT  FORMULATION 


A.  STIFFNESS  MATRIX 

1.  Shell  Element 

When  the  thickness  of  a  structure  is  small  compared  to  its  lateral  dimension, 
computational  efficiency  can  be  gained  by  using  a  two-dimensional  shell  element  in  the 
formulation.  A  nine-noded,  quadratic,  isoparametric  shell  element  [Ref  1]  was  selected 
for  use  in  the  finite  element  code  and  is  shown  in  Fig.  2.1.  Transverse  shear  deformation 
was  included  in  the  formulation  which  results  in  the  phenomenon  that  straight  lines 
originally  normal  to  the  undeformed  shell  midplane  surface  remain  straight,  but  not 
necessarily  normal  to  the  deformed  midplane  surface.  Deflections  were  assumed  to  be 
linear  functions  of  the  local  coordinate  in  the  thickness  direction.  One  advantage  of  this 
nine-noded  element  over  other  four  and  eight-noded  shell  elements  is  that  the  central 
nodal  point  is  very  convenient  when  selecting  a  representative  point  for  velocities  and 
pressures  on  the  surface  of  the  element  for  use  with  the  fluid-structure  interaction 


Figure  2.1  Nine  Noded  Isoparametric  Shell  Element  (From  Ref.  1) 
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Following  the  formulation  outlined  by  Bathe  [Ref.  1],  the  Cartesian  coordinates  of 
a  point  within  the  interior  of  the  undeformed  shell  element  with  local  coordinates  r,  s  and 
t  are  given  by 

°x(r,s,t)=l[0xl+|ak0V„l,]Nk(r,s) 

9  (  t  > 

0y(r,s,t)  =  S  °yk+-ak°Vnky  Nk(r,s)  (2.1) 

k=i  V  Z  J 

“z(r.s,t)=|;(0zk+iak0V„‘]Nl(r,s) 

where  °xk,  °yk  and  °zk  are  the  undeformed  Cartesian  coordinates  of  nodal  point  k,  ak  is  the 
shell  thickness  at  nodal  point  k,  °Vnk  is  the  unit  vector  normal  to  the  undeformed  shell 
mid-surface  at  nodal  point  k,  and  Nk  are  the  two-dimensional,  bi-quadratic  shape 
functions.  Similarly,  the  Cartesian  coordinates  of  a  point  within  the  deformed  shell  are 

9  f  t  \ 

1x(r,s,t)  =  XI 'xk+-ak'Vnkx jNk(r,s) 

k=l  v  ^  ' 

,y(r,s,t)  =  tf,y,+7aklvONk(r,s)  (2.2) 

k=l  v  y 

'z(r,s,t)  =  iflzk+tak'V,k,jNk(r;S) 
k=i  V  z  / 

where  'xk,  'yk  and  'zk  are  the  deformed  Cartesian  coordinates  of  nodal  point  k  and  1  Vk,  is 
the  unit  vector  normal  to  the  deformed  shell  mid-surface  at  nodal  point  k. 

The  shell  displacement  components  are  found  by  subtracting  the  undeformed 
coordinates  of  Eq.  (2.1)  from  the  deformed  coordinates  of  Eq.  (2.2)  to  yield 


where  V k = 1  Vk -° Vk  is  the  increment  in  the  direction  cosines  of 0  Vk . 
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The  two-dimensional  shape  functions  are  obtained  by  appropriately  multiplying 
together  the  one-dimensional,  isoparametric  shape  functions  which  correspond  to  the 
local  coordinates  r  and  s.  Assume  that  the  one-dimensional  elements  have  their  first  nodal 
point  located  at  r  =  1,  the  second  at  r  =  -1,  and  the  third  at  r  =  0.  The  quadratic  shape 
functions  for  these  three-noded  elements  in  the  local  r  coordinate  are 

N,r  =i(r2  +r)  ,  N„={(r!-r),  N„  =  l-r!  (2.4) 

Similarly,  the  one-dimensional  shape  functions  in  the  local  s  coordinate  are 

N^tfj’+s).  N2s  =  ;(s2  -s)  ,  N„  =  l-s!  (2.5) 

Note  that  th.e  k111  shape  function  is  unity  at  the  k'h  node  and  zero  at  all  the  other  nodes. 
Multiplying  these  shape  functions  together  yields  the  two  dimensional  shape  functions 
shown  in  Table  2. 1  below  where  the  local  coordinates  of  the  k1*1  node  have  been  included 
for  completeness. 


Nodal  Point 

r 

S 

Shape  Function 

1 

i 

l 

N,=NlrNls 

2 

-i 

l 

N2=N2rNls 

3 

-i 

-l 

N3  =  N2rN2s 

4 

i 

-l 

N4  =  NlrN2s 

5 

0 

l 

N5=N3rNls 

6 

-1 

0 

N6  =  N2rN3s 

7 

0 

-1 

N7=N3rN2s 

8 

1 

0 

N8=NlrN3s 

9 

0 

0 

N9  =  N3rN3s 

Table  2.1  Shell  Element  Shape  Functions  (see  Fig.  2.1) 


A  convenient  way  to  represent  the  direction  cosines  increment,  Vk  is  to  define 
two  unit  vectors, 0  V,k  and  °"\£,  that  are  orthogonal  both  to  0  Vk  and  to  each  other  as  shown 
in  Fig.  2.1 .  Let  Vk  be  a  linear  combination  of  °Vk  and  °V£  such  that  Vk  =  °Vk(3k-0V2kak 
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where  ak  and  pk  represent  the  small  angle  rotations  of  the  undeformed  normal  vector 0  Vnk 
about  the  vectors  °V,k  and  °\£,  respectively.  Since  their  is  no  unique  set  of  °Vkand  °Y* 
orthogonal  to  °Vk,  define  °Vk  =  ey  x  °Vk  and  let  °\£  =  °Vk  x  °Vk  where  ey  is  a  unit 
vector  parallel  to  the  y-axis.  For  the  case  that  °Vk  is  parallel  to  ey,  let  °Vk=e2.  Note  that 
the  unknown  direction  cosines  increment,  Vk,  has  now  been  replaced  with  two  unknown 
nodal  rotations  which  are  physically  much  more  meaningful  and  make  the  application  of 
slope  boundary  conditions  simpler.  Substituting  for  Vk  into  Eq.  (2.3)  yields 

u(r,s,t)  =  l(ul+ta1("Vlk<pi-«V‘,al))Nk(r,s) 
v(r,M)  =  i(v„  +^at("V1kPk-Xa1,)}jk(r,s)  (2.6) 

w(r,s,t)  =  X(wi1  +^ak(°V1kp|t-0V2kzak)jNk(r,s) 

In  order  to  calculate  the  element  strains,  partial  derivatives  of  the  displacements 

5u 

need  to  be  determined.  For  example;  to  find  the  longitudinal  strain  sxx  =  — ,  the  partial 

ox 

derivatives  of  u  in  the  local  (r,s,t)  coordinates  are  undertaken  first,  then  they  are 
premultiplied  by  the  inverse  of  the  Jacobian  matrix,  J,  to  determine  the  partial  derivatives 
with  respect  to  the  Cartesian  (x,y,z)  coordinates.  The  displacement  derivatives  in  local 
coordinates  are  given  by 


r  > 
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0 
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5w 
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<  wk  • 

5w 

0 

0  0 
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<*k 

Jdt. 

_ 

k. 

where  gk  =  -~ak°  V2k  and  gk=kak0V,k.  Note  that  the  matricies  in  Eq.  (2.7)  have 

dimension  3  x  45  for  the  nine-noded  shell  element. 

The  Jacobian  matrix  is  given  by 


9x 

5y 

5z 

Sr 

5r 

5r 

9x 

5y 

5z 

9s 

5s 

5s 

9x 

5y 

5z 

.aT 

5t 

5t. 

(2.8) 


where  the  components  in  J  are  determined  by  taking  partial  derivatives  of  Eq.  (2.1) 
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The  partial  derivatives  of  the  shape  functions  in  Eq.  2.9  can  easily  be  found  as  products  of 
the  one-dimensional  shape  functions  and  their  derivatives  in  both  the  r  and  s  local 
coordinates.  The  results  are  shown  in  Table  2.2  with  the  one-dimensional  shape  function 
derivatives  determined  by 


dNtr 

dr 

dNls 


ds 


r  +  t, 


s  +  f, 


3N2r 

5r 

dN2s 


ds 


r 


1 

2  5 


dN3 

dr 

dN 


3s 


ds 


=  -2r 

=  -2s 


(2.10) 


k 

dNk 

dNk 

dr 

ds 

1 

dN,r 

N 

dr 

N  *'■ 
lr  ds 

2 

dN2r 

2  N 
dr  ls 

^2r  ds 

3 

dN2r 

2  N 
dr  N2s 

*2r  ds 

4 

dNlr 
dr  "2s 

N,  ^ 

,r  ds 

5 

dr  N" 

N 

a 

6 

a  Nj- 

N 

ds 

7 

dr  u 

N 

”  a 

8 

3N"n 

-  a  * 

N  **- 

"  a 

9 

dr  "3s 

N  ** 

3r  ds 

Table  2.2  Shell  Element  Shape  Function  Derivatives 


The  displacement  derivatives  in  Cartesian  coordinates  for  a  point  in  the  shell  element 
specified  with  local  (r,s,t)  coordinates  then  become 
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(2.10) 
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du 
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ds 

dv 

<5w 

dw 

dx 

du 

dy 

du 

’  =  UV1' 
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ds 
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>  < 

9 

dx 

dv 

dy  i 

dv 

¥ 

dv 

ds 

dv 

. 

dx 

dw 

dy 

dw 

MJ]-1' 

~dr 

dw 

"ds" 

dw 

^dz> 

.  dz- 

.  dt , 

.  dz  - 

L  dt  . 

which  allows  the  strain-displacement  matrix  B  to  be  directly  assembled  where  B  has 
dimension  6  x  45  for  the  nine-noded  shell  element  and  is  defined  below. 


(2.11) 


The  objective  now  is  to  determine  the  constitutive  matrix  relating  the  element 
strains  to  stresses,  a  =  Dshel)  e.  If  the  shell  material  is  made  up  of  orthotropic,  thin 
laminae,  then  there  are  three  coordinate  systems  to  be  concerned  with:  the  global  (x,y,z), 
local  (r,s,t)  and  material  (1,2,3)  as  shown  in  Fig.  2.2.  Here  it  is  assumed  that  the  principal 
material  axis,  usually  the  one  with  the  highest  elastic  modulus,  is  aligned  in  the  “1” 
direction  and  the  local  “t”  axis  is  coincident  with  the  material  “3”  axis.  The  global 
coordinate  axes  are  typically  aligned  along  a  major  dimension  of  the  structure. 


X 

Figure  2.2  Coordinate  Systems  for  a  Shell  Element 


The  constitutive  relation  for  an  orthotropic  lamina  aligned  in  the  material  (1,2,3) 
coordinate  system  is  given  by  a,  =  Dortho  s,  which  in  expanded  form  becomes  (Ref.  2) 

Qn  Q 12  0  0  0  0  £n 

^22  Q 12  Q22  0  0  0  0  S22 

a„  0  0  0  0  0  0  e33 

1  r=  n  r  (2.12) 

x,2  0  00  q66  0  0  y  ,2 

-c23  0  0  0  0  kQ44  0  y23 


where  the  effects  of  shear  deformation  in  the  shell  have  been  included.  The  normal  stress 
and  strain  in  the  thickness  direction,  a3  and  e3,  are  neglected  in  this  shell  formulation.  The 
Qy  factors  are  related  to  the  material  constants  by 

^  Eu  ^  U2|E|1  rv _ E22 _  /->  1  “5\ 


E « 1  u,  |  E .  * 

1  u12u21  1  Ul2U2I 


Q22  - 


1  u12o21 


(2.13) 


Q44  -  G23  , 


Q55  -  EI31  > 


Q66  -  Gl 


10 


where  E  is  the  modulus  of  elasticity,  v  the  Poisson  ratio,  and  G  the  shear  modulus.  The 
shear  correction  factor,  k  =  f ,  accounts  for  the  fact  that  the  finite  element  representation 
of  the  shell  yields  a  constant  transverse  shear  stress,  whereas  it  actually  has  a  non¬ 
constant  distribution  when  calculated  with  continuum  mechanics. 

The  next  step  is  to  transform  the  constitutive  law  from  the  material  coordinate 
system  to  the  local  coordinate  system.  The  material  constitutive  matrix  Dortho  can  be 
converted  into  the  local  coordinate  constitutive  matrix  D,ocal  with  a  tensor  transformation 


of  the  form 


<*r  =  Dlocal  £r  =  T'1  Dortho  T  Sr 

Dlnoa|  =  T  Dortho  T 


'local 


(2.14) 


where  ar  =  [an  ctss  att  xrs  xst 
matrix  T  is  given  by  [Ref.  1] 


er  =  [Srr  Sss  %  Yrs  Yst  Yrtl  and  the  transformation 


T  = 


1? 

w2 

mi 

n? 

l,m, 

l,n, 

122 

__  2 

in  2 

n2 

l2m2 

m2n2 

l2n2 

^2 

m3 

_2 

n3 

l3m3 

m3n3 

l3n3 

21,  h 

2m,m2 

2n,n2 

l,m2  +  l2m, 

mtn2  -hm2n1 

l,n2  +  l2n, 

21*2 1 3 

2m2m3 

2n2n3 

l2m3  +l3m2 

m2n3  +  m3n2 

l2n3+l3n2 

21, 13 

2m,m3 

2n,n3 

l,m3  +l3m, 

m{n3  +  m3nl 

l,n3  +l3n, 

1, 

=  cos(eF , 

e,), 

m,  =  cos(ej 

»®i)»  ni 

=  cos(et,e,) 

12 

II 

O 

O 

C/> 

✓ - s 

fb 

>®2)> 

m2  =  cos(e§ 

)®2  ^2 

=  cos(et,e2) 

13 

=  cos(er- , 

■  «3)> 

m3  =  cos(e5 

,e3),  n3 

=  cos(et,e3) 

(2.15) 


where 


The  unit  vectors  eF  and  e5  have  been  “orthogonalized”  according  to 

ec  xe, 


(2.16) 


er  = 


|es  x  e, 


e,  =  e,  x  ef 


(2.17) 


since  er,  es,  and  et  are  in  general  not  orthogonal.  Assume  that  the  material  coordinate 
system  is  rotated  counter-clockwise  at  an  angle  0  with  respect  to  the  local  coordinate 
system  (with  the  material  3  axis  coincident  with  the  local  t  axis).  Let  m=cos(0)  and 
n=sin(0).  The  direction  cosines  of  Eq.  (2.16)  reduce  to  li=m,  I2-  -n,  mi=n,  m2=m, 
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I3=m3=n,=n2=n3=0.  Substituting  T  from  Eq.  (2.15)  and  Dortho  from  Eq.  (2.12)  into  Eq. 
(2.14)  yields 

Qn  Q12  o  q16  o  o 

Q|2  Q22  ®  Q26  0  0 

0  0  0  0  0  0 

D" =  Q,«  Q»  0  q«  0  _o 

0  0  0  0  kQ44  Q45 

0  0  0  0  Q45  kQ55_ 

where 

Qi,  -Q, im4+2(Q12+2Q66)mV+Q22n4 
Q22  =  Qun4+2(Q12+2Q66)m2n2+Q22m4 

Q44  =  Q44m2+Q55n2 

Q55  =  Q44n2+Q55m2 

Q66  =  (Qll+Q22‘2Ql2'2Q66)m2n2+Q66(m  +n  ) 

Qn  =  (Qii+Q22-4Q66)m2n2+Q12(m4+n4) 

Qi6  =  (Q 11  -Qi2-2Q66)m3n+(Qi2-Q22+2Q66)mn 

Q26  =  (QirQi2‘2Q66)mn3+(Qi2‘Q22+'2Q66)m  n 

Q45  =  (Qs5-Q44)mn 

The  next  step  is  to  transform  the  constitutive  matrix  from  local  coordinates  to 
global  (x,y,z)  coordinates.  This  is  done  with  another  tensor  transformation  of  the  form 

cfx  —  ® shell  -  T  DIocai  T  sx 

D«,  =  T1  Dm  T  (2.20) 

T  T 

where  crx  =  [cr^  (Tyy  xXy  XyZ  xX2]  and  sx—  [exx  Syy  yXy  YyZ  yxz]  and  the  tensor 
transformation  matrix  T  is  also  given  by  Eq.  (2.15).  The  direction  cosines  for  this  case  are 
determined  by 


(2.18) 


(2.19) 
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1,  =  cos(ex,eF), 

m,  =  cos(ey,eF), 

n,  =  cos(ez,eF) 

12  =cos(ex,e,). 

m2  =  cos(ey,e5), 

n2  =cos(ez,e3)  (2.21) 

13  =cos(ex,et), 

m3  =  cos(ey,et), 

n3  =  cos(ez,et) 

We  finally  have  enough  information  now  to  generate  the  element  stiffness  matrix 

This  matrix  relates  element  displacements  and  rotations  to  the  applied  element  forces  and 
moments  according  to 

kelem  <*  =  kem  (2-22) 

where  d  is  the  vector  of  nodal  displacements  and  rotations  and  has  45  components.  fdem 
is  the  vector  of  applied  nodal  forces  and  moments  and  also  has  45  components.  ke,em  is 
determined  by  the  volume  integral 

k«,em  =  jBTDshellBdV  (2.23) 

V 

where  dV  =  |J|drdsdt.  The  integration  is  carried  out  numerically  using  the  Gauss 
quadrature  method.  Three  integration  points  are  needed  in  the  r  and  s  coordinates,  while 
two  points  are  sufficient  in  the  t  coordinate. 

2.  Solid  Element 

When  components  of  a  structure  are  not  thin  enough  to  justify  use  of  the  shell 
element  formulated  above,  a  three-dimensional  solid  element  must  be  used  in  the  finite 
element  analysis.  The  27  node,  quadratic,  isoparametric,  solid  (brick)  element  was  chosen 
for  use  in  the  present  finite  element  code  and  is  shown  in  Fig.  2.3.  Athough  it  is  a 


Figure  2.3  27-Noded  Isoparametric  Solid  Element 
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computationally  expensive  element  compared  to  other  fewer  noded  solid  elements,  it  has 
the  advantage  of  providing  easy  compatibility  when  mixed  with  the  9-noded  shell 
element. 

The  development  of  the  solid  element  stiffness  matrix  closely  parallels  the 
development  of  the  stiffness  matrix  for  the  shell  element  mentioned  previously.  Consider 
the  element  of  Fig.  2.3  in  the  local  coordinate  system.  The  Cartesian  coordinates  of  a 
point  within  the  interior  of  the  solid  element  are  given  by 

x(r,s,t)  =  X  *kNk,  y(r,s,t)  =  X  YkNk’  z(r,s,t)  =  X  zkNk  (2-24) 

k=l  k=l  k=l 

where  xk,  yk  and  zk  are  the  undeformed  Cartesian  coordinates  of  nodal  point  k  and  Nk  are 
the  three-dimensional,  quadratic  shape  functions.  Similarly,  the  displacement  components 
of  a  point  within  the  solid  element  are  given  by 

27  V  27 

u(r,s,t)  =  Z  ukNk,  v(r,s,t)  =  £  vkNk,  w(r,s,t)  =  £  wkNk  (2.25) 

k=l  k=l  k=l 

where  uk,  vk  and  wk  are  the  displacements  of  nodal  point  k  (three  displacement  DOF  per 
node).  The  three-dimensional  shape  functions  are  obtained  by  appropriately  multiplying 
together  the.  one-dimensional  shape  functions  as  was  done  previously  for  the  shell 
element.  The  resulting  shape  functions  are  shown  in  Table  2.3  along  with  the  local 
coordinates  of  the  nodal  points.  The  one-dimensional  shape  functions  in  the  r  and  s 
coordinate  systems  were  defined  in  Eqs.  (2.4)  and  (2.5)  while  in  the  t  coordinate  system 

N„=i(t2+t)>  N2t  =  y(t2  - 1)  and  N3t=l-t2  (2.26) 

Table  2.4  summarizes  the  shape  function  partial  derivatives.  The  derivatives  in  the  r  and 
s  coordinates  were  previously  given  by  Eq.  2.10  while  the  derivatives  in  the  t  coordinate 
are  given  by 


5Nu 

dr 


dN2r 
dr  F 


and 


5N3T 

dr 


=  -2r 


(2.27) 
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Nodal  Point 

r 

s 

t 

Shape  Function 

1 

i 

■■ 

MM 

N,=NtrNlsNlt 

2 

-i 

19 

19 

N2=N2rNlsNIt 

3 

-i 

-1 

19 

N3  =  N2rN2sNlt 

4 

i 

-1 

19 

N4  =  NlrN2sNlt 

5 

0 

1 

1 

N5=N3rNlsNlt 

6 

-1 

0 

1 

N6=N2rN3sNlt 

'  7 

0 

-1 

1 

N7  =  N3rN2sNlt 

8 

1 

0 

1 

N8  =  NlrN3sNlt 

9 

0 

0 

1 

N9  =  N3rN3sNlt 

10 

1 

1 

2 

N10  =  NlrNlsN2t 

11 

-1 

1 

2 

N„  =N2rN|SN2t 

12 

-1 

-1 

2 

N12  =  N2rN2sN2t 

13 

1 

-1 

2 

N'13  =  NlrN2sN2t 

14 

0 

1 

2 

N14  =  N3rNlsN2t 

15 

-1 

0 

2 

N15  =  N2rN3sN2t 

16 

0 

-1 

2 

Ni6“N3rN2sN2t 

17 

1 

0 

2 

Nlv  =  NlrN3sN2t 

18 

0 

0 

2 

NI8  =  N3rN3sN2t 

19 

1 

1 

3 

N,9  =  NlrNlsN3t 

20 

-1 

1 

3 

N20=N2rNlsN3t 

21 

-1 

-1 

3 

N21=N2rN2sN3t 

22 

1 

-1 

3 

N22  =  NlrN2sN3t 

23 

0 

1 

3 

N23  =  N3rNlsN3t 

24 

-1 

0 

3 

N24  =  N2rN3sN3t 

25 

0 

sip 

3 

N25  =  N3rN2sN3t 

26 

1 

3 

N26  =  NlrN3sN3t 

27 

0 

99 

3 

N27  =  N3rN3sN3t 

Table  2.3  Solid  Element  Shape  Functions 
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In  order  to  generate  the  strain  displacement  matrix,  B,  we  first  need  to  find  the 
displacement  derivatives  in  local  coordinates  and  then  perform  a  transformation  to 
Cartesian  coordinates.  This  transformation  is  accomplished  through  use  of  the  Jacobian 

matrix.  The  elements  in  the  Jacobian  matrix,  J,  defined  in  Eq.  (2.8)  can  be  found  by 
taking  partial  derivatives  of  Eq.  (2.17)  according  to 


dx  ^  dNk 
dr  ~hX'  dr  ’ 

ay  f  aNk 
dr  k=1  k  ar  ’ 

dz  ^  aNk 
dr~hZk  dr 

dx  y^  SNk 
ds  ~  i“?Xk  ds  ’ 

ay  f  aNk 
as  £rk  as  ’ 

dz  ^  aNk 
as  k=1  Zk  as 

(2.28) 

5x  ^  3Nk 

a'trXk  at  ’ 

dy  y1  dNk 
dt  ~£tyk  dt  ’ 

Qz  9N  k 

di~hZk  at 

The  local  coordinate  displacement  derivatives  are  determined  by  differentiating 
the  displacements  given  in  Eq.  2.25.  The  results  are  shown  in  matrix  form  as 


au 

aNk  A  A 

0  0 

,  -v 

ai7 

ar 

Uk 

au 

~ds 

ii 

lM" 

aNk  n  A 
— k  0  0 
as 

K 

<  v  r 

v  k 

au 

aNk  ^  A 

Wk, 

— 

— L  0  0 

,  at . 

.  at  J 

dv 

'  aNk  “ 

a7 

ar 

Uv 

dv 

ds 

■-£ 

k=l 

o 

'i  ® 
© 

k 

vk  * 

dv 

5Nk 

0  0 

Wk. 

M. 

dt 

aw 

aNk' 

0  0  — — 

r  " 

ar 

dr 

Uk 

aw 
*  ~a7 

27 

-2 

k=l 

o 

o 

K 

'vk 

aw 

n  n  aNk 

Uk. 

0  0  ~z~~ 

,  dt 

dt  . 

(2.29) 
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Premultiplying  the  matricies  in  Eq.  (2.29)  by  the  inverse  of  the  Jacobian  matrix  gives  the 
displacement  derivatives  in  Cartesian  coordinates 


(2.30) 


da 

da 

'dv' 

r  > 

dv 

dw 

dw 

dx 

dr 

dx 

a^ 

dx 

dT 

da 

da 

dv 

dv 

dw 

dw 

1  — 

i  “7 

>  < 

% 

— 

8s 

dy 

ds 

dy 

ds 

8a 

da 

dv 

dv 

dw 

dw 

.a. 

.dz. 

.  at . 

.  dz . 

,  at  > 

which  are  used  to  assemble  the  size  6x81  strain-displacement  matrix,  B,  defined  by 


yy 


du 

r  > 

dx 

V, 

dv 

w. 

dy 

t 

u2 

dw 

N 

‘ 

>  1 
fo  1 

au 

•  =  [B]« 

v2 

W2 

dx 

dy 

u , 

dw 

dv 

3 

dy  + 

dz 
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V27 

.  dx  + 
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lYxz 


The  next  task  is  to  find  the  stress-strain  matrix,  Dbrick.  For  the  solid  element  we 
assume  the  material  is  isotropic.  Thus,  no  tensor  transformation  is  necessary  to  convert 
the  constitutive  matrix  from  local  to  Cartesian  coordinates.  The  stress-strain  matrix  is 


(2.31) 


D 


E(l-o) 


brick 


(1  +  UX1-2U) 


1 

o 


u 


o 


1-0 

o 


l-o  1-0 
o 

1 


0 

0 


o 


1-0  l-o 
0  0 


l-o 
1  0 
l-2o 


0  0 

0  0  0  0 

where  E  is  the  modulus  of  elasticity  and  v  is  the  Poisson  ratio. 


0 

0 
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2(l-o) 

0 


0 

0 

0 

0 

l-2o 

2(l-o) 

0 


0 

0 

0 

0 

0 

l-2o 
2(1 -o)J 


(2.32) 
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The  shell  element  stiffness  matrix  can  now  be  found  with  the  volume  integral 

keien,  =  f  B T  D  brick  B  d  V  (2.33) 

V 

where  dV  =  |j|drdsdt.  The  integration  is  carried  out  numerically  using  the  Gauss 
quadrature  method.  Three  integration  points  are  needed  for  each  of  the  r,  s  and  t 
coordinates. 

3.  Assembly  into  Global  Stiffness  Matrix 

Once  the  stiffness  matrix  is  generated  for  an  element,  the  next  step  is  to  assign  the 
components  of  that  matrix  into  the  global  stiffness  matrix,  Kj.  The  relationship  between 
the  position  of  an  entry  in  the  element  matrix  and  the  corresponding  position  in  the  global 
matrix  is  obtained  by  conducting  a  table  look-up  in  a  connectivity  table  which  is 
generated  separately. 

Due  to  symmetry  and  the  highly  banded  nature  of  K^,  an  efficient  way  to  store  the 
components  of  it  is  through  use  of  the  skyline  storage  method  [Refs.  1  &  9].  This  scheme 
stores  in  a  large  vector  only  those  entries  which  are  within  the  “skyline”  of  the  upper  half 
of  Ks.  The  heights  of  the  columns  within  the  “skyline”  are  found  by  examining  the 
connectivity  table  for  the  global  DOF  of  interest  and  then  locating  the  element  attached  to 
this  global  DOF  which  has  the  largest  spread  of  global  DOF  associated  with  it.  The 
height  of  the  column  is  the  difference  between  the  largest  and  smallest  DOF  associated 
with  this  element. 

B.  MASS  MATRIX 

1.  Shell  Element 

In  order  to  simplify  the  computations  of  generating  the  shell  element  mass  matrix, 
a  lumped  mass  model  was  used.  This  method  entails  assigning  a  fraction  of  the  total 
element  mass  to  each  element  node.  This  model  has  the  added  benefit  of  producing  a 
diagonal  mass  matrix  which  tremendously  reduces  the  storage  requirements  and  the 
numerical  operations  involved  when  solving  for  the  structural  dynamic  response. 
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The  matrix  O,  whose  diagonal  elements  consist  of  the  fractions  of  the  total 
element  mass  carried  by  each  node,  is  found  by  integrating  over  the  area  of  the  shell 
element  in  local  coordinates  according  to 

(D=pajNNTdA  (2.34) 

A 

where  N  is  the  column  vector  of  shape  functions,  dA  =  |j|dr  ds ,  a  is  the  shell  thickness 
and  p  is  the  mass  density,  assumed  constant  throughout  the  element.  If  the  integration  is 
carried  out  using  the  three-point  Lobatto  rule  where  the  integration  points  correspond  to 
the  nodal  points  [Ref.  9],  then  the  off-diagonal  terms  of  <f>  become  zero.  For  the  case  of  a 
rectangular  shell  with  constant  thickness,  the  Jacobian  matrix  is  constant  throughout  the 
shell.  This  produces  diagonal  components  of  O  equal  to 

t  ^-[l  1  1  1  4  4  4  4  16] 

where  Ve  is  the  total  volume  of  the  element.  Thus,  the  central  node  receives  j  of  the  total 
shell  mass,  each  edge  node  receives  i  of  the  total  and  each  comer  node  of  the  total  for 
the  case  of  a  rectangular  element. 

The  mass  contribution  of  each  shell  node  is  then  assigned  to  the  diagonal 
locations  in  the  element  mass  matrix  corresponding  to  each  of  the  three  translational 
degree-of-freedoms  (DOF)  for  that  node.  The  remaining  two  diagonal  locations  in  the 
matrix  for  each  node  correspond  to  moment  of  inertia  terms  which  are  neglected  in  this 
formulation  (they  are  assigned  small,  non-zero  values  in  order  to  keep  the  mass  matrix 
non-singular).  The  result  is  a  diagonal  element  mass  matrix  which  can  be  stored  as  a 
vector  of  length  45. 

2.  Solid  Element 

The  lumped  mass  model  was  also  used  to  determine  the  element  mass  matrix  for 
the  solid  element.  The  matrix  <1>,  whose  diagonal  elements  consist  of  the  fraction  of  the 
total  element  mass  carried  by  each  node,  is  found  by  integrating  over  the  volume  of  the 
solid  element  in  local  coordinates  according  to 


20 


(2.35) 


<D=  p  Jn  Nt  dV 

V 

where  dV  =  |j|drdsdt  and  p  is  the  mass  density,  assumed  constant  throughout  the 
element.  For  the  case  of  a  parallelepiped  shaped  element,  the  Jacobian  matrix  is  constant 
throughout  the  element.  This  yields  diagonal  components  of  O  equal  to 

£^\\  l  l  l  4  4  4  4  16  1  1  1  1  4  4  4  4  16  4  4  4  4  16  16  16  16  64l 
216 1 

where  Ve  is  the  total  volume  of  the  element.  Thus,  the  node  at  the  center  of  the  element 
receives  jf  of  the  total  element  mass,  each  face  node  receives  jf  of  the  total,  each  edge 
node  receives  k  of  the  total  and  each  comer  node  -jk  of  the  total  element  mass.  The  mass 
contribution  of  each  node  is  then  assigned  to  the  diagonal  locations  in  the  element  mass 
matrix  corresponding  to  each  of  the  three  translational  DOF  for  that  node.  The  result  is  a 
diagonal  element  mass  matrix  which  can  be  stored  as  a  vector  of  length  8 1 . 

3.  Assembly  into  Global  Mass  Matrix 

Since  the  element  mass  matricies  are  already  diagonalized,  their  assembly  into  the 
global  mass  matrix,  Ms,  is  easily  performed.  The  components  of  the  element  mass 
matricies  are  directly  assigned  into  positions  in  the  global  mass  matrix  as  specified  by  the 
connectivity  table.  The  use  of  a  diagonal  global  mass  matrix  requires  very  little  storage 
and  provides  great  computational  efficiency  in  the  transient  response  calculations. 

C.  DAMPING  MATRIX 

The  system  damping  matrix  is  formed  as  a  linear  combination  of  the  global  mass 
and  stiffness  matricies,  termed  Rayleigh  damping,  as  follows 

Cs=aMs  +  (5Ks  (236) 

where  a  and  p  are  constants.  If  a=0  and  P  is  nonzero  then  the  damping  tends  to  filter  out 
high  frequency  components.  Setting  P=0  and  a  nonzero  filters  out  low  frequency 
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components.  When  a  and  P  are  both  nonzero  then  the  damping  is  minimum  at 
intermediate  frequencies  and  becomes  large  at  both  low  and  high  frequencies. 
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III.  FLUID-STRUCTURE  INTERACTION 


A.  GOVERNING  EQUATIONS 

When  a  structure  is  submerged  in  a  fluid  and  subjected  to  an  incident  acoustic 
wave,  the  resulting  structural  motion  affects  the  strength  of  the  scattered  pressure  wave. 
In  addition,  changes  in  the  scattered  pressure  strength  affect  the  dynamic  response  of  the 
structure.  Hence,  the  structural  motion  is  coupled  to  the  scattered  fluid  pressure.  The 
governing  equations  for  the  fluid-structure  response  can  thus  be  represented  as  a  set  of 
coupled,  linear,  ordinary  differential  equations. 

Application  of  Newton’s  second  law  to  a  structure  submerged  in  an  infinite 
acoustic  medium  results  in 

M,x  +  C,j  +  K,x  =  -CAf(p,+ps)  (3.1) 

where  Ms,  Cs  and  Ks  are  the  structural  mass,  damping  and  stiffness  matrices, 
respectively,  x  is  the  structural  displacement  vector,  Af  is  a  diagonal  matrix  of  the  fluid 
element  areas,  pj  is  the  vector  of  incident  pressures  and  ps  is  the  vector  of  scattered 
pressures.  G  is  the  pressure/force  compatibility  matrix  which  assigns  the  fluid  pressures 
acting  on  the  wet  structure  fluid  nodes  to  forces  acting  on  the  structural  nodes.  The  G 
matrix  also  relates  the  fluid  particle  velocities  and  structural  nodal  velocities  according  to 

Gtx  =  u,+us  (3.2) 

where  U[  is  the  incident  normal  fluid  particle  velocity  and  us  is  the  scattered  (reflected) 
normal  fluid  particle  velocity. 

The  scattered  pressure  in  Eq.  (3.1)  can  be  obtained  by  applying  a  doubly  asymptotic 
approximation  (DAA)  which  is  based  upon  the  representation  of  surface  motion  as  a 
linear  combination  of  orthogonal  fluid  boundary  modes.  The  simplest  approximation, 
DAAb  as  given  by  Geers  [Ref.  1 8]  is 

Mfps  +  pcAfps  =  pcMfus  (3.3) 

where  Mf  is  the  fluid  mass  matrix,  p  is  the  fluid  density  and  c  is  the  speed  of  sound  in 
the  fluid.  Solving  Eq.  (3.2)  for  us  and  substituting  into  Eq.  (3.3)  yields 

Mfp,  +  pcAfps  =  pcMf(GTx-u,)  (3.4) 
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Eqs.  (3.1)  and  (3.4)  are  the  governing  equations  for  the  fluid-structure  interaction 
problem  with  unknown  variables  x  and  ps. 

B.  FLUID  MASS  MATRIX  DEVELOPMENT 

The  fluid  mass  matrix,  Mf ,  is  produced  by  using  a  boundary  element  formulation 
to  solve  Laplace’s  equation  for  an  incompressible,  inviscid,  irrotational  and  infinite  fluid 
subjected  to  motion  of  the  structure’s  wet  surface.  Mf  is  determined  solely  from  the 
geometrical  characteristics  of  the  finite-element  representation  of  the  wet  surface  of  the 
body  and  is  fully  populated.  The  governing  equations  for  fluid  motion  are  [Ref  4] 

V2<p  =  0  in  the  fluid  domain 

—  =  -u  on  the  structure' s  surface  (3.5) 

5n 

Vcp  =  0  infinitely  far  from  the  structure 

where  <p  is  the  velocity  potential,  n  is  a  unit  vector  normal  to  the  structure  surface  and  u  is 
the  normal  velocity  of  the  structure’s  surface. 

A  simple  source  formulation  [Ref.  4]  can  be  used  to  solve  Eq.  (3.5)  according  to 

<P(P)  =  l““TdS(q) 

s  P’q  (3.6) 

f(p)=-2MP)+jo(q);;eq(p’q)ds(q) 

on  s  r  (p,q) 

where  r(p,q)  is  the  magnitude  of  the  vector  r(p,q)  from  the  receiver,  point  p,  to  the 
transmitter,  point  q;  a(q)  is  the  source  strength;  S  is  the  wetted  structural  surface;  and  0 
is  the  angle  between  r(p,q)  and  n(p).  The  geometry  is  shown  in  Fig.  3.1.  The  kinetic 


Figure  3.1.  Source  Formulation  Geometry 
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energy  of  the  fluid  can  be  expressed  as 

Tf  =-ypJ^(p)9(p)dS(p)  (3.7) 

S  ^ 


where  p  is  the  fluid  density. 

If  it  is  assumed  that  the  source  strength  is  constant  throughout  each  fluid  element, 
then  Eqs.  (3.6)  and  (3.7)  can  be  written  in  matrix  form  as 


<P(P)  =  B(p)o 

(3.8) 

!^(P)  =  -C(p)a 

(3.9) 

Tf=-2pj^(p)  dS(p)cp(p) 

(3-10) 

where  the  elements  of  B  and  C  are  given  by 

h/  ,  r dS(qj) 

Sj  rCP i  j  ) 

(3.11) 

rcos0(pi,qi) 

C||(P|) = 2,t5ii  *  s  "7(p3q^dSlqi) 

(3.12) 

Using  the  I]  scheme  described  in  Ref.  4,  Eqs.  3.11  and  3.12  can  be  simplified  to 


equal  to  (p').  Eliminating  a  in  Eq.  (3.8)  and  the  transpose  of  Eq.  (3.9)  yields 

on 
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(3.15) 


<p(p)  =  B(p)C-,(p')u 

l^(p)jt =-uTc’T(p,)CT<p) 

We  can  now  obtain  a  more  useful  expression  for  the  fluid  kinetic  energy  by 
substituting  Eq.  (3.15)  into  Eq.  (3.10)  to  produce 

T,  =|puTC-T(p'|jcI(p)dS(p)B(p)  C-'(p')u  (3.16) 

If  we  use  one  point  integration  for  the  integral  in  Eq.  (3.16)  with  p  =  p’  then  the  fluid 
kinetic  energy  can  be  expressed  as 

Tf  =|puTAfBC'u  (3.17) 

where  Af  is  a  diagonal  matrix  of  the  fluid  element  areas.  Let  E  =  AfBC_l  and  define  the 
fluid  mass  matrix  by 

Tf=|uTMfu  (3.18) 

Then  the  symmetric  fluid  mass  matrix  can  now  be  expressed  as 

Mf={p(E  +  ET)  (3.19) 


C.  STABILITY  CONSIDERATIONS 

The  solution  of  the  interaction  Eqs.  (3.1)  and  (3.4)  can  be  further  stabilized  with 
respect  to  the  choice  of  time  increment  by  injecting  one  of  the  coupled  equations  into  the 
other  as  recommended  by  DeRuntz  [Ref.  3].  Solving  Eq.  (3.1)  for  x ,  substituting  into  Eq. 
(3.4)  and  premultiplying  both  sides  of  the  resulting  equation  by  AfMf 1  yields 
Afps  +(/xAfMf,Af)ps  =  -  pcAfGTM;'Csx-  pc AfGTM;Ksx 
-(/xAfGTM;'GAf)Pl  -(/xAfGTM;'GAf)ps  -  pcAfu, 

Defining  Dn=/xAfMj1Af,  D  =/xAfGTM;‘GAf,  D!=Dn+Ds  and  D2=-pcAfGTM;' 
yields  the  augmented,  second  fluid-structure  interaction  equation 

AfPs  +  D,ps  =  D2(Csx  +  Ksx)-Dsp,-pcAfu,  (3.21) 

The  presence  of  the  u,  term  in  Eq.  (3.21)  presents  numerical  difficulties  when  an 
incident  shock  wave  is  modeled.  This  can  be  alleviated  for  a  spherical  incident  wave  by 
defining  a  modified  pressure  vector 
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Pm~Ps  +  ^Pi  (3.22) 

where  T  is  a  diagonal  matrix  of  the  direction  cosine  elements  y,  and  y^cos  (<t>j).  Oj  is  the 
angle  between  a  vector  coincident  with  R,,  the  distance  from  the  origin  of  the  incident 
spherical  wave  to  the  ith  fluid  node  on  the  wet  surface  of  the  structure,  and  nt,  the  wet 
surface  normal  vector  of  the  ith  fluid  node  as  shown  in  Fig.  3.2. 


STRUCTURE 

Figure  3.2.  Charge  Geometry 


The  wet  surface  pressure,  pTot,  is  now  given  by 

Prot  =Pi+Ps  =  (l-r)p,+pM  (3-23) 

Solving  Eq.  (3.22)  for  ps  and  substituting  into  Eq.  (3.21)  produces  the  modified, 
augmented  fluid-structure  interaction  equations 

Msx  +  Csx+Ksx  =  -GAf(pM  +(l-r)Pl)  (3.24) 

AfpM+D,pM=  D2(Csx  +  Ksx)+(D,r-Ds)pI+Afrp,-pcAfuI  (3.25) 
For  a  spherical  wave  in  a  homogeneous  fluid,  the  incident  wave  pressure  at  the  ith  fluid 
node,  p,  ,  is  related  to  the  incident  wave  fluid  particle  velocity  at  the  i*  fluid  node,U[. , 
and  the  incident  wave  pressure  at  the  standoff  distance,  p1?  by  the  following 

*,,(t>=4^(,>+i-p''(,))  (3-26) 

S  (  R;  -  S') 

p'(t)“iTp'lt-— J  (327) 
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where  R  is  a  diagonal  matrix  of  the  radii  R(  and  S,  the  standoff  distance,  is  defined  to  be 
the  distance  from  the  origin  of  the  spherical  wave  to  the  nearest  point  on  the  wet 
structure.  Converting  Eq.  (3.26)  into  vector  format  and  rearranging  to  solve  for  p,  yields 

p,  =  pcr''u, -cR'lp1  (3.28) 

which,  when  substituted  into  Eq.  (3.25),  causes  the  terms  involving  ii,  to  cancel  out 
producing 

AfpM+D,pM=  D2(Csx  +  Ksx)  +  (D,r-Ds-cAfR-!r)p1  (3.29) 

Defining  H  =  D,T-DS -cAfR''r  yields 

AfpM+®iPM=  D2(Csx  +  Ksx)-Hpj  (3.30) 

D.  NUMERICAL  SOLUTION  OF  THE  FLUID-STRUCTURE 
INTERACTION  EQUATIONS 

The  solution  of  the  modified,  augmented,  fluid-structure  interaction  Eqs.  (3.24)  and 
(3.30)  is  obtained  by  using  a  staggered  solution  procedure  in  conjunction  with  the  central 
finite  difference  method.  Assuming  the  structure  is  initially  at  rest  with  initial  conditions 

-At 

x°  =  x°  =  x~  2  =  Pm  =  0  ^  assuming  the  incident  pressure,  pI?  is  known  at  all  times, 
the  first  calculation  is  to  find  the  initial  acceleration  vector.  Solving  Eq.  (3.24)  for  x  at 
time  zero,  then  integrating  for  velocity  and  displacement  for  the  first  time  step  yields 

x°  =  -M''(GAf(p„  +(I-r)pf)+Csx°  +  Ksx°)  =  -M;'GAt(I-r)p! 

At  At 

X  2  =  X  2  +X%t  =  X°At 

XA*  =X0+X^At  =  X^At  (3.31) 

Applying  Eq.  (3.30)  over  the  first  time  step  yields 

At  At  (  At  At  \  At 

ArpJ  +D,pJ  =  D!|C,xt+Ksi>  J-Hp,=  (3.32) 

Substituting  =~(pm  -Pm)  and  pi  =  |(p„ +P„')  into  Eq.  (3.32)  produces 
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TAf(PM -Pm)  +  JDi(p°m+Pm)  =  D2 


At 


f  At  AtA  At 

Csx 2  +  Ksx 2  -Hp,2 

V  J 


(3.33) 


which,  after  solving  for  the  modified  pressure  at  time  At,  gives 

(] 

VAt 


Pm  = 


\  Wi  a  (\  yi  f  At  ^  (j  yi  * 

-Af+jD,]  ^Af-iD,yM4^Af+]D1J  D2^+K5x^J-^Af+|D,j  Bp}  (3.34) 


1  1 


At 


-1 


1  1 


Let  D3  =  —  Af+—  D,  I  5  D4-  .  Af—  D,,  D5-D3D4,  D6-D3D2,  D7— -D3H 


At 


(  At  aA 


then 


p^D^+D^  +Ksx^ 


+  D~p,2 


(3.35) 


The  displacements  and  pressures  for  subsequent  time  steps  are  calculated  in  a  similar 
manner  with  the  difference  equations  for  the  nth  time  step  reducing  to 

xnAt  =-M;‘(GAf(p^  +(l-r)p^))  +  CsxnAt  +Ksxn 
(•♦i)4*  (-1)“ 


_  nAt 


=  X  +XnA,At 


■+-  At 

X("+l,At  =  XnA*  +X  At 


(3.36) 


(•4)“ 

2 

(  (  1 


X  =  —  (x"“'  +  X 


nAt  ,  _(n+l)At 


) 


p5T  =DjP«  +d» 


V 


(**0*  R)a 

Csx  +KjX 


+D7p, 
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IV.  CODE  APPLICATIONS 


A.  VERIFICATION  EXAMPLES 

In  order  to  verify  the  accuracy  of  the  finite  element  code,  several  example 
problems  involving  a  wide  variety  of  structural  geometries  and  loadings  are 
demonstrated.  The  exact  analytical  solutions  are  also  presented  for  comparison  when 
available. 

1.  Free  Vibration  of  a  Dry,  Isotropic  Rod 

Consider  the  transient  response  of  an  isotropic  rod  subjected  only  to  an  initial  load 
which  is  released  at  time  zero.  One  end  of  the  rod  is  rigidly  fixed  while  the  other  end  is 
free  as  shown  in  Fig.  4.1.  The  finite  element  model  was  constructed  with  a  uniform  mesh 
of  20  brick  elements  laid  end-to-end.  For  the  case  of  a  steel  rod  with  length  L=10  ft, 
cross-sectional  area  A=1  ft2,  Young’s  modulus  E=4.32xl09  lb/ft2,  density  p=15.2  slug/ft3 
and  force  F=1000  lbs,  the  transient  response  is  shown  in  Fig.  4.2. 


z 


The  analytical  solution  for  the  free  response  of  a  rod  with  a  constant  cross- 
sectional  area  and  Young’s  modulus  can  be  determined  by  the  wave  equation: 
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(4.1) 


52u  2  d2u 

at7"'  a7 


where  u  is  the  axial  displacement  at  a  point  in  the  rod. 


is  the  speed  of  sound  in 


the  rod  material.  The  solution  to  Eq.  4.1  is  obtained  by  using  the  method  of  separation  of 
variables  with  the  result  given  by  [Ref.  10] 


u(x,t) 


8FL  #-«  (-1)"  (2n+l)7ix  (2n  +  l)7tct 

— ; - -rsin — COS — 

7t2AE^(2n+l)2  2L  2L 


(4.2) 


Substituting-in  values  for  the  current  problem  reduces  Eq.  4.2  to 

00  (_l)2n 

u(L,t)  =  1.88x1 0"6  ft  Er - rrcos[(2n+  l)2648t]  (4.3) 

„=o  (2n  + 1) 

for  displacement  at  the  rod  tip.  The  analytical  solution  including  the  first  five  terms  in  the 
series  is  also  plotted  in  Fig.  4.2.  Excellent  agreement  occurs  between  the  FEM  and 
analytical  solutions. 


Time  (msec) 


Figure  4.2.  Dry  Rod  Transient  Response 
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2.  Dry,  Isotropic  Plate  Subjected  to  a  Step  Pressure  Loading 

Now  consider  the  dynamic  response  of  a  flat  plate  subjected  to  transverse  loading 
such  that  plate  bending  occurs.  A  simply  supported,  isotropic,  square  plate  subjected  to  a 
uniformly  distributed  step  loading  was  modeled.  For  simplicity,  only  one-quarter  of  the 
plate  was  modeled  with  two  sides  simply  supported  and  symmetric  boundary  conditions 
enforced  on  the  remaining  sides,  as  shown  in  Fig.  4.3.  A  5  x  5  element  mesh  composed  of 
nine-noded  shell  elements  was  used  in  the  finite  element  solution.  The  transient  response 
at  the  midpoint  of  a  square  aluminum  plate  with  length  a=5  ft,  thickness  t=T  in., 
E=1.483xl09  lb/ft2  and  p=5.25  slug/ft3  subjected  to  a  1  lb/ft2  uniformly  distributed  step 
loading  is  shown  in  Fig.  4.4.  Sinusoidal  motion  occurs  with  an  amplitude  of  6.6xl0‘5  ft. 
and  a  period  of  18.6  msec. 


v=0 


Figure  4.3.  Dry  Plate  Problem  Geometry 

The  analytical  solution  for  the  plate  transient  response  can  be  obtained  by 
examining  both  the  plate  static  response  and  the  undamped  natural  frequencies  of 
vibration.  For  a  simply-supported  square  plate  subjected  to  a  uniform  load  P0,  the  static 
deflection  at  the  midpoint  is  given  by  [Ref.  1 1] 


33 


/  \  sm  1 

16P0  ^  ^  (-1)  2 

^  _r  6  pv  /hvJ  2^  .  >2  /  \2l^ 

n  D  m=1'3n=l-3rnn[(7)  +(7) 

(4.4) 

Et3 

where  D  -  ,  ,  x 

12(l-u2) 

(4.5) 

This  series  converges  rapidly.  If  only  the  first  four  terms  in  the  series  are  retained  and 
v=0.3,  then  Eq.  4.4  can  be  simplified  to 

P0a4  ^ 

w  =  0.0443  3  (4.6) 

Substituting  in  values  for  the  current  problem  yields  the  static  deflection  at  the  plate 
midpoint  of  3.3xl0'5  ft.  For  a  plate  subjected  to  a  uniform  step  loading,  the  plate  midpoint 
deflection  is  expected  to  oscillate  about  the  static  deflection  value  with  a  maximum 
amplitude  equal  to  twice  the  static  deflection.  This  is  exactly  what  occurs  in  the  FEM 
solution. 


’1 
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Figure  4.4.  Dry  Plate  Transient  Response 
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The  undamped  natural  frequencies  of  an  isotropic,  simply  supported  square  plate 
are  given  by  [Ref.  12] 

(m2  +  n2)  (4.7) 

For  a  uniformly  distributed  step  loading,  the  first  mode  of  vibration  is  expected  to  be  the 
dominant  mode.  Thus,  setting  m=n=l  and  substituting  the  constants  for  the  current 
problem  into  Eq.  4.7  yields  coj [=334.6  rad/sec.  This  corresponds  to  a  vibration  period  of 
18.8  msec  which  is  commensurate  with  the  results  shown  in  Fig.  4.4. 

3.  Wet,  Isotropic  Rod  Subjected  to  a  Step  Pressure  Loading 


®m,n  = 


71  D 


In  order  to  assess  the  dynamic  response  of  a  submerged  structure  to  an  acoustic 
shock  wave,  consider  a  rod  fixed  at  one  end  and  subjected  to  a  fluid  pressure  wave  at  the 
other  end  as  shown  in  Fig.  4.5.  Symmetric  boundary  conditions  are  imposed  on  the  sides 
of  the  rod  to  produce  a  one-dimensional  response  along  the  longitudinal  axis  only. 

z 


unit 


t 

pressure  wave 


Figure  4.5.  Wet  Rod  Problem  Geometry 


If  the  structure  is  homogeneous  with  a  constant  cross-sectional  area,  the  incident 
pressure  wave  will  propagate  through  the  rod  at  the  speed  of  sound  for  that  material,  c0. 
This  speed  is  related  to  the  material  density  and  elastic  modulus  by 
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(4.8) 


When  the  pressure  wave  reaches  the  fixed  boundary,  it  will  be  completely  reflected  as 
another  compression  wave  of  the  same  magnitude,  causing  the  resultant  stress  at  points 
immediately  behind  this  reflected  wave  to  be  double  the  value  of  stress  created  by  the 
original  incident  wave  [Ref.  13].  Upon  reaching  the  free  edge,  the  wave  will  again  be 
reflected.  However,  since  the  characteristic  impedance,  pc,  of  the  fluid  is  much  less  than 
that  of  the  structure,  this  reflected  wave  will  be  a  tensile  wave.  If  the  surrounding  fluid  is 
air,  whose  impedance  is  negligible  compared  to  that  of  the  structure,  then  points  in  the 
structure  immediately  behind  this  reflected  tensile  wave  will  have  a  stress  value  equal  to 
the  incident  wave  pressure.  When  the  tensile  wave  reaches  the  fixed  boundary,  it  will  be 
completely  reflected  as  another  tensile  wave  of  the  same  magnitude.  Thus  points 
immediately  behind  this  reflected  tensile  wave  will  have  a  zero-stress  state.  The 
theoretical  stress  response  for  a  point  at  the  midpoint  of  a  20  ft  long  aluminum  rod  in  air 
subjected  to  a  1  lb/ft2  incident  wave  is  shown  in  Fig.  4.6. 


0  4  8  12  16 
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Figure  4.6.  Theoretical  Rod  Midpoint  Stress  Response 
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If  the  rod  is  surrounded  with  water,  then  the  structural  dynamics  is  altered 
significantly.  The  impedance  of  water  is  3600  times  the  impedance  of  air  which  allows 
the  stress  wave  energy  to  be  more  readily  transmitted  into  the  fluid  at  the  free  end  instead 
of  being  nearly  completely  reflected  as  was  the  case  with  air.  This  causes  a  dampening 
effect  on  the  structural  response.  Figure  4.7  shows  the  finite  element  solution  of  the 
longitudinal  normal  stress  time  history  at  the  midpoint  of  a  20  ft.  long,  aluminum  rod 
impinged  by  a  1  lb/ft2  incident  pressure  wave  with  three  different  rod  widths  (square  rod 
cross-section). 


0  2  4  6  8  10  12  14  16 
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Figure  4.7.  Wet  Rod  Longitudinal  Stress  History 

The  aspect  ratio,  AR,  is  defined  as  the  rod  length  divided  by  rod  width.  As  the  cross- 
sectional  area  of  the  rod  becomes  larger  while  the  rod  length  is  held  constant,  the  fluid 
dampening  has  a  much  greater  effect  causing  the  stress  amplitude  to  decay  faster  over 
time.  In  addition,  the  added  fluid  mass  becomes  a  larger  fraction  of  the  total  system  mass 
as  the  aspect  ratio  becomes  smaller.  This  results  in  larger  periods  of  oscillation  with 
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smaller  aspect  ratios.  For  example,  if  the  cross-sectional  area  of  the  rod  is  doubled  while 
the  rod  length  is  held  constant,  then  both  the  structural  mass,  wetted  surface  area,  and 
structural  stiffness  in  the  longitudinal  direction  will  double.  However,  the  larger  wetted 
surface  can  now  entrain  a  deeper  depth  of  fluid  which  causes  the  virtual  fluid  mass  to 
more  than  double.  Thus  the  total  system  mass  grows  faster  than  the  stiffness  as  the  aspect 
ratio  is  lowered  causing  the  system  frequency  of  oscillation  to  drop. 

By  contrast,  the  transient  response  of  a  rod  in  a  vacuum  subjected  to  a  step  load  at 
the  free  end  shows  no  dependence  on  the  rod  width  as  shown  in  Fig.  4.8.  The  rod 
geometry  and  material  are  identical  to  the  wet  rod  examined  above.  In  this  case,  if  the  rod 
width  is  doubled  then  both  the  total  system  mass  and  longitudinal  stiffness  double  which 
produce  no  change  in  the  undamped  natural  frequency.  The  dynamic  response  shown  in 
Fig.  4.8  compares  favorably  to  the  theoretical  response  shown  in  Fig  4.6. 


Figure  4.8.  Dry  Rod  Longitudinal  Stress  History 
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4.  Wet,  Isotropic  Plate  Subjected  to  a  Step  Pressure  Loading 

If  the  plate  discussed  in  part  IV.A.2  is  submerged  in  water  and  subjected  to  a  unit 
pressure  wave,  the  displacement  of  the  plate  midpoint  behaves  as  shown  in  Fig.  4.9.  The 
maximum  displacement  is  bounded  by  almost  twice  the  static  deflection  shown  in  Fig. 
4.4,  but  the  response  here  is  very  highly  damped.  Furthermore,  the  period  of  oscillation 
for  the  wet  plate  is  roughly  three  times  the  period  of  the  dry  plate  due  to  the  presence  of 
the  added  fluid  mass.  As  a  matter  of  fact,  the  fluid  mass  contributes  more  to  the  total 
system  mass  than  does  the  structural  mass. 


Figure  4.9.  Wet  Plate  Midpoint  Transient  Response 

The  wet  plate  midpoint  pressure  component  time  histories  are  shown  in  Fig.  4.10. 
The  total  pressure  is  simply  the  sum  of  the  incident  and  scattered  pressure  and  is  the 
actual  pressure  felt  by  a  particle  immediately  behind  the  scattered  wave.  As  can  be  seen, 
the  total  pressure  oscillates  about  the  incident  wave  magnitude  with  a  period  of 
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oscillation  identical  to  that  of  the  plate  midpoint  displacement  response.  The  total 
pressure  never  drops  below  zero  here  because  water  cannot  support  a  tension  wave  the 
way  a  solid  material  can. 


Figure  4.10.  Wet  Plate  Fluid  Pressure  Response 

5.  Static  Deflection  of  an  Axisymmetric  Cylinder 

In  order  to  verify  the  accuracy  of  finite  element  code  when  curved  surfaces  are 
encountered,  consider  the  static  response  of  a  cylindrical  shell  clamped  at  one  end  and 
with  a  radial  distributed  load  at  the  free  end  as  shown  in  Fig.  4.11.  Only  one  quarter  of 
the  cylinder  was  modeled  due  to  symmetry  considerations.  The  finite  element 
formulation  used  240  shell  elements  in  a  10x24  mesh  with  a  greater  concentration  of 
elements  located  near  the  free  end  where  the  displacement  gradient  was  highest.  The 
radial  deflection  is  plotted  in  Fig.  4.12  along  with  the  exact  solution  [Ref.  14].  Excellent 
agreement  occurred  between  the  two  solutions. 
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Fixed  Boundary 


Figure  4.1 1.  Axisymmetric  Cylinder  Problem  Geometry 
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Figure  4.12.  Radial  Deflection  of  the  Axisymmetric  Cylinder 


6.  Static  Deflection  of  an  Angle-Ply  Composite 

To  check  the  accuracy  of  the  code  when  anisotropic  material  effects  are  included, 
consider  the  static  response  of  a  two-layer,  square,  antisymmetric  (+0,-0),  angle-ply 
composite  such  as  that  shown  in  Fig.  4.13. 


Figure  4.13.  Angle  Ply  Composite 

.  nx  .  Try  . 

The  composite  is  transversely  loaded  according  to  p  =  p0  sin — sin —  where  a  is  the 

3.  3 

edge  width.  The  material  properties  for  each  ply  are  as  follows: 

E,/E2=25,  G12/E2=0.6,  G22/E2=0.5,  v12=0.25 

When  the  finite  element  code  was  run  on  a  one-quarter  model  using  a  5x5  shell  element 
mesh  with  a=10  ft.,  E,=2.5062xl09  lb/ft2,  p0=100  lb/ft2 ,  h=a/6  and  various  ply  angles,  the 
normalized  midpoint  deflections  were  as  shown  in  Fig.  4.14.  The  published  solution  was 
also  obtained  numerically  [Ref.  15].  The  error  between  the  two  solutions  is  mainly 
attributed  to  lumping  the  load  distribution  over  each  element  at  the  central  element  node 
in  the  present  formulation.  This  force  lumping  method  was  used  because  the  same 
computer  code  was  used  to  solve  fluid-structure  interaction  problems  and  with  those  type 
of  problems,  only  one  element  nodal  point  is  available  to  apply  the  pressure  loading.  If  a 
truly  consistent  force  vector  was  used  for  this  dry,  composite  problem,  then  slightly  better 
results  would  be  expected. 
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Figure  4.14.  Midpoint  Deflection  for  a  Two  Layer  Angle-ply  Composite 


7.  Wet,  Isotropic  Sphere  Subjected  to  a  Step  Pressure  Loading 

The  last  verification  example  examines  the  dynamic  response  of  a  submerged, 
unrestrained,  spherical  shell  to  an  incident  plane  wave  as  shown  in  Fig.  4.15. 


Incident 
Plane  Wave 


y 


Figure  4.15.  Spherical  Shell  Problem  Geometry 
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The  finite  element  formulation  uses  an  upper-half  model  consisting  of  24  shell 
elements  arranged  in  a  hemisphere  with  symmetric  boundary  conditions  applied  along  the 
equator.  The  sphere  has  a  radius  a=1.0  m  and  thickness  b=0.02  m.  The  sphere  is  made  of 
steel  with  properties  E=206.84  GPa,  v=0.33,  p=7784.5  kg/m3  and  the  surrounding  water 
has  sound  speed  c=  146 1.2  m/s  and  density  p=999.6  kg/m3.  The  incident  wave  of  pressure 
1 .0  N/m2  is  modeled  as  a  spherical  wave  with  a  very  large  standoff  distance  in  order  to 
approximate  a  plane  wave.  The  normalized  velocity  response  of  the  points  on  the  sphere 
closest  and  farthest  from  the  charge  are  shown  Fig.  4.16  and  Fig.  4.17,  respectively.  The 
analytical  solution  was  obtained  using  the  the  method  of  separation  of  variables  [Ref.  19] 
and  is  shown  for  comparison.  A  numerical  solution  obtained  with  the  USA/DYNA3D 
software  combination  [Ref  20]  is  also  shown.  Since  the  USA/D YNA3D  software  uses  the 
more  exact  DAA2  doubly  assymptotic  approximation  which  includes  the  effects  of 
curvature  in  the  fluid  mass  matrix  construction,  it  is  expected  to  yield  slightly  better 
results  than  the  present  formulation  which  is  based  on  the  simpler  DAA1  method. 
However,  the  present  method  tracks  reasonably  well  with  the  analytic  solution  for  the 
early  time  response. 

/ 
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Normalized  Velocity  (p  cu/P,)  Normalized  Velocity  (p  cu/P,) 


Figure  4.17.  Far  Element  Response 


B.  UNBALANCED  SANDWICH  COMPOSITE  IMPACT 
STUDY 

Sandwich  composites  are  becoming  increasingly  more  popular  as  structural 
components  for  military  applications.  These  composites  are  constructed  of  two  strong 
facesheets  separated  by  a  thick,  lightweight,  weaker  core.  The  resulting  composite  can 
have  a  much  higher  specific  strength  than  a  comparable  monolithic  material.  One 
application  of  this  composite  of  interest  to  the  Navy  is  its  use  in  the  Advanced 
Performance  Mast  System  (APMS).  This  particular  sandwich  composite  has  titanium  and 
glass  fiber  reinforced  plastic  (GRP)  facesheets  separated  by  a  phenolic,  Nomex  fiber 
reinforced  honeycomb  core.  Since  the  facesheets  are  of  different  material,  this  type  of 
sandwich  composite  is  termed  “unbalanced.” 

The  objectives  of  this  study  are  to  generate  a  numerical  model  of  the  unbalanced 
sandwich  composite  and  then  use  this  model  to  simulate  the  dynamic  response  of  the 
composite  subject  to  impact  loading,  with  and  without  material  failure. 

An  experimental  impact  study  was  conducted  by  Fuller  [Ref.  7]  on  samples  of  the 
unbalanced  sandwich  composite  used  on  the  APMS.  The  samples  had  length  =  12.0  in., 
width  =  2.75  in.  and  a  nominal  thickness  of  1.18  in.  The  three-point  impact  tests  used  a 
fixture  which  held  the  composite  beams  in  a  simply  supported  configuration  with  the 
impactor  striking  the  center  of  the  beam  as  shown  in  Fig.  4.18.  The  supports  were  located 
0.5  in.  from  each  end.  Five  strain  gages  were  bonded  to  the  beam  as  shown  in  Fig.  4.19. 


Impactor  Force 


Figure  4.18.  Unbalanced  Sandwich  Composite 
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opposite  side 


F igure  4.19.  Strain  Gage  Placement 

Strain  gage  #5  was  located  at  the  center  of  the  beam  on  the  opposite  side  of  the  impact 
with  the  remaining  strain  gages  located  at  the  beam  quarter  points. 

The  first  test  to  be  modeled  is  the  impact  of  a  125  lb.  mass  dropped  from  a  height 
of  1.0  in.  (0.0254  m.)  onto  the  midpoint  of  the  sandwich  composite.  This  drop  height  was 
small  enough  that  material  failure  did  not  occur.  The  experimental  test  used  a  force 
transducer  located  between  the  beam  and  impactor  to  measure  the  impact  force 
throughout  the  impact  process.  Fig.  4.20  shows  the  impact  force  time  history  which  was 


Figure  4.20.  Experimental  Impact  Force  Time  History  for  1.0  in.  Drop 


47 


used  as  an  input  to  the  finite  element  code.  In  the  experimental  setup,  a  thin  strip  of  brass 
was  fastened  to  the  center  of  the  beam  to  spread  out  the  contact  load  over  the  width  of  the 
beam.  This  was  modeled  in  the  finite  element  formulation  by  treating  the  impact  load  as  a 
line  force  across  the  width  of  the  beam. 

The  finite  element  model  of  the  sandwich  composite  was  constructed  with  12 
nine-noded  shell  elements  of  equal  length  laid  end-to-end.  Each  element  consisted  of 
three  layers  stacked  on  top  of  each  other  with  different  material  properties  in  each  layer. 
Table  4.1  summarizes  the  material  properties  for  each  layer. 


En  (psi) 
E22  (psi) 
Gj2  (psi) 
G,3  (psi) 
G23  (psi) 


Titanium 


15.5x10 
15.5  xlO6 
5.77xl06 
5.77x1 06 
5.77xl06 
0.342 


Honeycomb 


p  (slug/ft3)  9.18 

Nominal  thickness  (in.)  0.1 

Actual  thickness  (in.)  0. 1 06 


Table  4.1.  Sandwich  Composite  Material  Properties 

The  experimental  strain  gage  readings  for  the  first  impact  test  are  shown  in  Fig. 
4.21  while  the  numerical  simulation  results  are  shown  in  Fig.  4.22.  Since  the  beam  is 
symmetric  and  no  failure  occurs,  the  strain  histories  for  strain  gages  #4  &  #6  are  identical 
as  are  the  strains  for  gages  #2  &  #3.  The  numerical  simulation  tracks  reasonably  well 
with  the  experimental  results  for  strain  gages  #2  &  #3  which  are  located  on  the  top  of  the 
composite  and  undergo  compression.  However,  the  numerical  simulation  differs  from  the 
experimental  results  for  the  strain  gages  located  on  the  beam  bottom.  If  we  assume  that 
the  beam  behaves  according  to  linear  elastic  beam  theory  and  neglect  the  shear 
deformation,  then  the  strain  at  the  midpoint  of  the  beam  (gage  #5)  should  always  be  a 
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Figure  4.21.  Experimental  Results  of  The  Impact  Test  with  Drop  Height  =  1.0  in. 
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Figure  4.22.  Numerical  Simulation  of  the  Impact  Test  with  Drop  Height  =  1 .0  in, 


factor  of  2.2  greater  than  the  strain  at  the  beam  quarterpoints  (#4  &  #6).  This  factor  of 
2.2  accounts  for  the  supports  not  being  located  exactly  at  the  beam  ends.  If  the  supports 
were  exactly  at  the  beam  ends  than  the  factor  would  be  2.0  as  expected.  The  numerical 
simulation,  which  is  based  on  linear  elastic  theory,  does  indeed  show  the  strain  at  the 
beam  midpoint  to  be  approximately  a  factor  of  2.2  higher  than  the  strains  at  the 
quarterpoints.  By  contrast,  the  experimental  test  shows  the  midpoint  strain  to  be  roughly  a 
factor  of  3.1  times  the  quarterpoint  strains.  This  disparity  cannot  be  accounted  for  by  the 
presence  of  the  shear  deformation  alone.  Thus,  although  the  numerical  simulation  can 
reasonably  predict  the  sandwich  composite  response  using  linear  elastic  theory,  some 
nonlinear  effects  are  clearly  present  in  the  actual  dynamics  of  the  sandwich  composite 
and  need  to  be  appropriately  modeled. 

The  next  experimental  test  to  be  simulated  uses  the  same  material  and  setup  as  the 
previous  case  but  with  an  impactor  drop  height  of  1.5  in.  This  height  is  large  enough  to 
cause  a  partial  failure  of  the  sandwich  composite  material.  The  experimentally 
determined  impact  force  time  history  is  shown  in  Fig.  4.23.  The  sharp  drop  in  the  impact 
force  at  time  0.01  seconds  is  due  to  the  composite  suddenly  relaxing  due  to  partial 
material  failure. 

The  experimental  strain  gage  readings  for  the  second  impact  test  are  shown  in  Fig. 
4.24  while  the  numerical  simulation  results  are  shown  in  Fig.  4.25.  The  strains  in  both  the 
experimental  test  and  simulation  up  to  the  time  of  failure  follow  the  same  trends  as  they 
did  for  the  impact  test  with  a  drop  height  of  1.0  in.  At  time  0.01  seconds,  the  sudden  drop 
in  the  impact  force  tends  to  cause  the  strain  at  all  points  in  the  beam  to  decrease.  This  can 
be  seen  in  Fig.  4.24  for  strain  gages  #3,  #4  &  #6.  However,  material  failure  occurs  at  the 
beam  quarterpoint  where  strain  gages  #2  &  #6  are  located.  This  causes  the  material  to 
relax  at  that  point  which  tends  to  cause  a  large  increase  in  tensile  strain  for  strain  gage  #6 
and  a  large  increase  in  compressive  strain  for  strain  gage  #2.  The  increase  in  strain  at  the 
quarterpoint  due  to  material  failure  is  much  more  significant  than  the  decrease  in  strain  at 
the  quarterpoint  due  to  the  lower  impactor  force.  The  end  result  is  that  the  strain  increases 
sharply  at  the  beam  quarterpoint  at  the  time  of  failure. 
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Figure  4.23.  Experimental  Impact  Force  Time  History  for  1.5  in.  Drop 

In  order  to  simulate  failure  in  the  finite  element  analysis  of  the  sandwich 
composite,  the  modulus  of  elasticity  of  each  layer  for  the  two  elements  bordering  the 
beam  quarter  point  was  reduced  by  40%  at  the  time  of  failure.  The  global  stiffness  matrix 
was  then  regenerated  and  the  transient  analysis  continued.  As  can  be  seen  in  Fig.  4.25,  the 
strain  response  can  be  reasonably  predicted  using  this  method.  An  alternate  method  to 
simulate  failure  would  be  to  reduce  the  effective  load  carrying  thickness  of  a  particular 
layer  where  failure  initiates.  Although  these  methods  may  seem  somewhat  arbitrary  in 
their  approach  to  fracture,  the  purpose  of  this  study  was  to  observe  the  effects  of  fracture 
on  the  dynamic  response  of  the  composite.  Further  research  needs  to  be  devoted  to 
determine  how  failure  initiating  in  the  core  material  weakens  the  structure  on  a  global 
scale  to  produce  the  results  seen  in  the  experimental  study. 
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Figure  4.24.  Experimental  Results  of  the  Impact  Test  with  Drop  Height  =  1 .5  in. 
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Figure  4.25.  Numerical  Simulation  of  the  Impact  Test  with  Drop  Height  =1.5  in. 


C.  EFFECTS  OF  COMPOSITE  LAYER  SMEARING  ON 
THE  FLUID-STRUCTURE  INTERACTION 

When  the  number  of  layers  in  a  composite  laminate  becomes  large,  the  computer 
time  required  to  generate  the  structural  stiffness  matrix  can  become  prohibitively 
expensive.  A  method  which  combines  the  material  properties  of  several  adjacent 
composite  laminae  into  an  single  layer  of  equivalent  stiffness,  termed  “smearing”,  is 
demonstrated  below  for  structures  excited  by  underwater  shock. 

Consider  a  filament  wound  S-glass/epoxy  composite  with  ten  equally  thick  layers 
layed-up  as  follows:  [0/4 5 /90/-4 5/90/4 5 /90/-45/90/90] .  This  fiber  orientation,  as  studied 
by  Rousseau  [Ref.  21],  is  intentionally  non-symmetrical  with  respect  to  the  midplane. 
The  conventional  method  to  generate  an  element  stiffness  matrix  of  this  composite  would 
be  to  first  determine  the  constitutive  relations  matrix,  Dortho,  for  each  layer  in  its 
associated  material  (1,2,3)  coordinate  system  (see  Fig.  2.2).  This  matrix  is  then 
transformed  into  global  coordinates  to  become  Dshell  and  used  along  with  the  strain- 
displacement  matrix,  Bshell,  to  determine  the  layer  stiffness.  The  stiffness  of  each  layer  is 
then  added  to  generate  the  element  stiffness  matrixfor  the  composite. 

The  method  used  here  to  reduce  the  computational  time  required  to  generate  the 
stiffness  matrix  involves  conducting  a  weighted  average  of  the  constitutive  matrix  for 
each  layer,  Dortho  layer  to  arrive  at  the  smeared  constitutive  matrix,  D0rth0i  smeared,  according 
to 

n 

^  i  ®  ortho ,  layeij 

®  ortho,  smeared 

i=l 

where  n  is  the  number  of  layers  and  hi  is  the  thickness  of  Ith  layer.  For  example,  if  it  is 
desired  to  smear  all  ten  layers  of  the  composite  mentioned  above  into  one  layer,  the 
smeared  constitutive  matrix  would  be  calculated  by 

D  ortho,  smeared  =  To(®0  +  ®90  +  2D4J  +2D_45)  (4.10) 
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This  smeared  constitutive  matrix  needs  to  be  calculated  only  once  for  the  structure, 
assuming  the  layer  thicknesses  and  material  characteristics  do  not  change  from  one 
element  to  the  next.  The  effective  material  constants  for  this  smeared  layer  can  also  be 
generated  from  the  smeared  constitutive  matrix,  if  desired  for  comparison  purposes,  but 
this  is  not  required  to  implement  the  smearing  formulation. 

To  illustrate  this  smearing  process  in  a  fluid-structure  interaction  problem, 
consider  a  simply  supported  plate,  5.0  ft.  on  each  side  and  is  1.0  in.  thick,  fabricated  with 
the  ten  layer  composite  material  discussed  above.  The  plate  is  exposed  to  water  on  one 
face  and  is  excited  by  a  1000  lb/ft  plane  incident  wave.  The  plate  midpoint  transverse 
deflection  time  history  is  shown  in  Fig.  4.26  and  the  midpoint  wet  surface  pressure  time 
history  is  shown  in  Fig.  4.27  for  various  smearing  scenarios. 

The  first  case  is  the  baseline  transient  response  with  all  composite  layers  used  to 
calculate  the  stiffness  matrix  in  the  conventional  manner.  The  next  case  smears  all  ten 
layers  into  one  effective  layer.  This  results  in  a  poor  approximation  to  the  unsmeared 
model  as  the  peak  midpoint  displacement  undershoots  by  1 1%  and  the  peak  displacement 
is  also  significantly  displaced  toward  an  earlier  time.  The  peak  surface  pressure  and  time 
to  peak  pressure  are  not  altered  significantly,  however.  The  third  case  uses  three  layers  in 
a  0  /  smear  /  90  orientation  where  the  smeared  layer  consists  of  the  eight  central  plies  of 
the  original  ten  layer  laminate.  Much  more  satisfactory  results  occur  with  this  case  as  the 
peak  midpoint  deflection  undershoots  by  only  3.5%  compared  to  the  unsmeared  model. 
In  addition,  the  midpoint  surface  pressure  response  is  nearly  identical  to  the  unsmeared 
case.  The  last  case  uses  five  layers  in  a  0  /  45/  smear  /  90  /  90  orientation  where  the 
smeared  layer  consists  of  the  six  central  plies  of  the  original  ten  layer.  Little  is  gained 
using  five  layers  instead  of  only  three  layers  as  the  peak  displacement  overshoots  by  only 
3.3%  and  the  surface  pressure  history  is  again  nearly  identical  to  the  unsmeared  model. 
Note  -  since  the  10  layer  composite  we  have  been  using  here  as  an  example  contains  two 
45  /  90/  -45  groups,  it  behaves  as  a  quasi-isotropic  material.  If  another  composite  layer 
orientation  had  been  selected  which  gave  a  much  greater  preference  to  one  direction  only, 
then  the  results  discussed  above  would  be  even  more  pronounced. 


54 


Composite  Plate  Midpoint  Deflection  Response 


Since  the  outer  layers  carry  most  of  the  load  in  thin  composite  structures  subjected 
to  transverse  loading  from  the  fluid-structure  interaction,  any  attempts  to  include  these 
layers  in  the  smearing  process  is  not  recommended.  On-the-other-hand,  composite  layers 
near  the  midplane  contribute  little  to  the  bending  stiffness  so  they  are  the  best  candidates 
for  smearing,  as  was  shown  in  Figs.  4.26  and  4.27. 
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V.  CONCLUSION 


The  general  purpose  finite  element/boundary  element  code  developed  here  has 
been  shown  to  be  useful  in  simulating  fluid-structure  interaction  and  impact  problems. 
Although  only  linear  elastic  effects  have  been  incorporated  at  this  point,  it  can  be  readily 
tailored  to  solve  specific  problems  such  as  those  dealing  with  modem  composite 
materials.  With  further  refinement,  it  will  prove  to  be  a  valuable  tool  in  future  underwater 
shock  research. 

The  following  items  are  recommended  to  be  incorporated  into  future  revisions  of 
the  code  to  improve  its  accuracy  and  scope: 

1.  Include  higher  order  methods  to  solve  the  integrals  of  Eqs.  3.11  and  3.12  in  the 
generation  of  the  fluid  mass  matrix  as  discussed  in  Ref.  4.  These  more 
advanced  methods  provide  better  results  with  coarse  meshes  and  handle  local 
curvature  effects  more  accurately. 

2.  Utilize  the  second  Doubly  Asymptotic  Approximation,  DAA2  [Ref.  3],  in  the 
governing  equations  for  the  fluid-structure  interaction.  This  method  is  a 
generalization  of  the  DAAj  to  a  symmetric  second  order  differential  equation 
with  improved  accuracy  in  the  intermediate  frequency  range.  It  can  better 
handle  global  curvature  of  the  submerged  structure. 

3.  Implement  the  Kwon  micromechanical  model  [Refs.  22-24]  to  determine  the 
micro-level  fiber  and  matrix  stresses.  These  stresses  can  be  used  in  a  stiffness 
reduction  correlations  to  better  estimate  material  property  degredation  and  to 
model  damage  progression  until  material  failure. 


59 


LIST  OF  REFERENCES 


1.  Bathe,  K.,  Finite  Element  Procedures  on  Engineering  Analysis ,  Prentice-Hall,  Inc., 
Englewood  Cliffs,  pp.  144-145,  pp.  194-297,  1982. 

2.  Vinson,  J.  R.  and  Chou,  T.,  Composite  Materials  and  their  Use  in  Structures,  John 
Wiley  &  Sons,  Inc.,  New  York,  pp.  201-223,  1975. 

3.  DeRuntz,  J.  A.,  “The  Underwater  Shock  Analysis  Code  and  its  Applications,”  60™ 
Shock  and  Vibration  Symposium  Proceedings,  Vol.  I,  pp.  89-107,  David  Taylor 
Research  Center,  Nov.,  1989. 

4.  DeRuntz,  J.  A.  and  Geers,  T.  L.,  “Added  Mass  Computation  by  the  Boundary  Integral 
Method,”  International  Journal  for  Numerical  Methods  in  Engineering,  Vol  12, 
1978,  pp.  531-550. 

5.  Huang,  H.,  Everstine,  C.  C.  and  Wang,  Y.  F.,  “Retarded  Potential  for  Analysis  of 
Submerged  Structues  Impinged  by  Weak  Shock  Waves,  ”  Computational  Methods  for 
Fluid-Structure  Interaction  Problems,  AMD- Vol  26,  American  Society  of 
Mechanical  Engineers,  December  1977. 

6.  Brasek.  T.  P.,  (1994)  Response  of  Dual-Layered  Structures  Subjected  to  Shock 
Pressure  Wave,  Master’s  Thesis,  Naval  Postgraduate  School,  Monterey,  California. 

7.  Fuller,  L.  B.,  (1994)  Damage  and  Compressive  Failure  of  Unbalanced  Sandwich 
Composite  Panels  Subject  to  a  Low-Velocity  Impact,  Master’s  Thesis,  Naval 
Postgraduate  School,  Monterey,  California. 

8.  Chiong,  L.  N.,  (1994)  Study  of  Air-Backed  and  Water-Backed  Composite  Cylinders 
Subjected  to  Underwater  Explosion,  Master’s  Thesis,  Naval  Postgraduate  School, 
Monterey,  California. 

9.  Akin,  J.  E.,  Finite  Elements  for  Analysis  and  Design,  Academic  Press,  pp.  359-421, 
1994. 

10.  Haberman,  R.,  Elementary  Applied  Partial  Differential  Equations,  2nd  Edition, 
Prentice-Hall,  Inc.,  pp.  113-124, 1987. 

11.  Ugural,  A.  C.  and  Fenster,  S.  K.,  Advanced  Strength  and  Applied  Elasticity,  2nd 
Edition,  Prentice-Hall,  Englewood  Cliffs,  pp.  419-421, 1987. 

12.  Timoshenko,  S.,  Theory  of  Plates  and  Shells,  McGraw  Hill,  New  York,  1959. 


61 


13.  Kolsky,  H.,  Stress  Waves  in  Solids,  New  York,  Dover  Publications,  Inc.,  1963. 

14.  Zienkiewicz,  O.  C.,  The  Finite  Element  Method,  3rd  Ed.,  McGraw-Hill,  London, 
p.361,  1977. 

15.  Whitney,  J.  M.  and  Pagano,  N.  J.,  “Shear  Deformation  in  Heterogeneous  Anisotropic 
Plates,”  Journal  of  Applied  Mechanics,  1970,  pp.  1031-1036. 

16.  Press,  W.  H.,  Teukolsky,  S.  A.,  Vetterling,  W.  T.,  Flannery,  B.  P.,  Numerical  Recipes 
in  FORTRAN.  The  Art  of  Scientific  Computing,  Cambridge  University  Press,  pp.  1- 
42,  1992. 

17.  Akai,  T.  J.,  Applied  Numerical  Methods  for  Engineers,  John  Wiley  &  Sons,  Inc., 
New  York,  pp.  289-313, 1994. 

18.  Geers,  T.  L.,  “Residual  Potential  and  Approximate  Methods  for  Three-Dimensionsl 
Fluid-Structure  Interaction  Problems,”  Journal  of  the  Acoustic  Society  of  America, 
Vol.  49, 1971,  pp.  1505-1510. 

19.  Huang,  H.,  “Transient  Interaction  of  Plane  Waves  with  a  Spherical  Elastic  Shell,” 
Journal  of  the  Acoustic  Society  of  America,  Vol.  45,  pp.  661-670,  1969. 

20.  Fox,  P.  K.,  (1992)  Nonlinear  Dynamic  Response  of  Cylindrical  Shells  Subjected  to 
Underwater  Side-On  Explosions,  Master’s  Thesis,  Naval  Postgraduate  School, 
Monterey,  California. 

21.  Rousseau,  M.  P.,  (1993)  Dynamic  Response  of  a  Filiment  Wound  Composite 
Cylinder  Exposed  to  Underwater  Shock,  Master’s  Thesis,  Naval  Postgraduate 
School,  Monterey,  California. 

22.  Kwon,  Y.  W.,  “Thermo-Elastoviscoplastic  Finite  Element  Plate  Bending  Analysis  of 
Composites,”  Engineering  Computations,  Vol.  9,  pp.  595-607,  July  1992. 

23.  Kwon,  Y.  W.,  “Calculation  of  Effective  Moduli  of  Fibrous  Composites  with 
Micromechanical  Damage,”  Composite  Structures,  Vol.  25,  pp.  187-192,  1993. 

24.  Berner,  J.  M.,  (1993)  Finite  Element  Analysis  of  Damage  in  Fibrous  Composites 
Using  a  Micromechanical  Model,  Master’s  Thesis,  Naval  Postgraduate  School, 
Monterey,  California. 


62 


INITIAL  DISTRIBUTION  LIST 


1.  Defense  Technical  Information  Center 
Cameron  Station 

Alexandria  VA  22304-6145 

2.  Library,  Code  052 
Naval  Postgraduate  School 
Monterey  CA  93943-5101 

3.  Professor  Young  W.  Kwon,  Code  ME/Kw 
Department  of  Mechanical  Engineering 
Naval  Postgraduate  School 

Monterey  CA  93943-5000 

4.  Professor  Young  S.  Shin,  Code  ME/Sg 
Department  of  Mechanical  Engineering 
Naval  Postgraduate  School 
Monterey  CA  93943-5000 

5.  Department  Chairman,  Code  ME/Kk 
Department  of  Mechanical  Engineering 
Naval  Postgraduate  School 
Monterey  CA  93943-5000 

6.  Naval  Engineering  Curricular  Office,  Code  34 
Naval  Postgraduate  School 

Monterey  CA  93943-5000 

7.  Erik  A.  Rasmussen 
Code  1720.4 

Naval  Surface  Warfare  Center,  Carderock  Division 
Bethesda  MD  20084-5000 

8.  Lt.  Robert  Burgio 
2913  Ordway  Drive 
Ellicott  City  MD  21042 


